#include<pmmintrin.h>

#include<time.h>

#include<xmmintrin.h>

#include<iostream>

#include <arm_neon.h>

#define N 64

using namespace std;



float** normal_gauss(float **matrix) //串行的高斯消去法

{

    for (int k = 0; k < N; k++)

    {

        float tmp =matrix[k][k];

        for (int j = k; j < N; j++)

        {

            matrix[k][j] = matrix[k][j] / tmp;

        }

        for (int i = k + 1; i < N; i++)

        {

            float tmp2 = matrix[i][k];

            for (int j = k + 1; j < N; j++)

            {

                matrix[i][j] = matrix[i][j] - tmp2 * matrix[k][j];

            }

            matrix[i][k] = 0;

        }

    }

    return matrix;

}



void SSE_gauss(float **matrix) //加了SSE并行的高斯消去法，4路

{

    __m128 t1, t2, t3, t4;

    for (int k = 0; k < N; k++)

    {

        float tmp[4] = { matrix[k][k], matrix[k][k], matrix[k][k], matrix[k][k] };

        t1 = _mm_loadu_ps(tmp);

        for (int j = N - 4; j >=k; j -= 4) //从后向前每次取四个

        {

            t2 = _mm_loadu_ps(matrix[k] + j);

            t3 = _mm_div_ps(t2, t1);//除法

            _mm_storeu_ps(matrix[k] + j, t3);

        }

        if (k % 4 != (N % 4)) //处理不能被4整除的元素

        {

            for (int j = k; j % 4 != ( N% 4); j++)

            {

                matrix[k][j] = matrix[k][j] / tmp[0];

            }

        }

        for (int j = (N % 4) - 1; j>= 0; j--)

        {

            matrix[k][j] = matrix[k][j] / tmp[0];

        }

        for (int i = k + 1; i < N; i++)

        {

            float tmp[4] = { matrix[i][k], matrix[i][k], matrix[i][k], matrix[i][k] };

            t1 = _mm_loadu_ps(tmp);

            for (int j = N - 4; j >k;j -= 4)

            {

                t2 = _mm_loadu_ps(matrix[i] + j);

                t3 = _mm_loadu_ps(matrix[k] + j);

                t4 = _mm_sub_ps(t2,_mm_mul_ps(t1, t3)); //减法

                _mm_storeu_ps(matrix[i] + j, t4);

            }

            for (int j = k + 1; j % 4 !=(N % 4); j++)

            {

                matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];

            }

            matrix[i][k] = 0;

        }

    }

}

void AVX(float **a){
    for(int k=0;k<N;k++){
        // float val[4] ={a[k][k],a[k][k],a[k][k],a[k][k]};
        // __m128 divsor =_mm_loadu_ps(val);
        __m256 divsor =_mm256_set1_ps(a[k][k]);
        int j;
        for (j = k + 1; j +8 < N; j+=8) {
            __m256 v =_mm256_loadu_ps(a[k]+j);
            v=_mm256_div_ps(v,divsor);
            _mm256_storeu_ps(a[k]+j,v);
        }
        for(;j<N;j++){
            a[k][j] = a[k][j] / a[k][k];
        }
        a[k][k]=1;
        for (izzt i = k + 1; i < N; i++) {
            // float val[4] = {a[i][k],a[i][k],a[i][k],a[i][k]};
            __m256 vaik=_mm256_set1_ps(a[i][k]);  
            int j;  
            for (j = k + 1; j +8< N; j+=8) {
                __m256 vakj=_mm256_loadu_ps(a[k]+j);
                __m256 vaij=_mm256_loadu_ps(a[i]+j);
                __m256 vx=_mm256_mul_ps(vakj,vaik);
                vaij =vaij -vx;
                // a[i][j] = a[i][j] - a[i][k] * a[k][j];
                _mm256_storeu_ps(a[i]+j, vaij);
            }
            for(;j<N;j++){
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            }
            a[i][k] = 0;
        }
    }

}






void Neon_gauss(float **matrix)
{
    float32x4_t t0,t1,t2,t3;
    for(int k=0;k<N;k++)
    {
        //float temp1[4]={a[k][k],a[k][k],a[k][k],a[k][k]};
        //t0=_mm_loadu_ps(temp1);
        t0=vld1q_dup_f32(matrix[k]+k);
        int j;

        for(j=k+1;j+3<N;j+=4)
        {
            t1=vld1q_f32(matrix[k]+j);
            t2=vdivq_f32(t1,t0);
            vst1q_f32(matrix[k]+j,t2);
        }
        for(;j<N;j++)
            matrix[k][j]/=matrix[k][k];

        matrix[k][k]=1.0;
        for(int i=k+1;i<N;i++)
        {
            //float temp2[4]={a[i][k],a[i][k],a[i][k],a[i][k]};
            t0=vld1q_dup_f32(matrix[i]+k);
            int j;
            for(j=k+1;j+3<N;j+=4)
            {
                t1=vld1q_f32(matrix[k]+j);
                t2=vld1q_f32(matrix[i]+j);
                t3=vmulq_f32(t0,t1);
                t2=vsubq_f32(t2,t3);
                vst1q_f32(matrix[i]+j,t2);
            }
            for(;j<N;j++)
                 matrix[i][j]-=matrix[i][k]*matrix[k][j];
            matrix[i][k]=0;
        }
    }

}

void print(float **matrix) //输出

{

    for (int i = 0; i < N; i++)

    {

        for (int j = 0; j < N; j++)

        {

            cout << matrix[i][j]<<" ";

        }

        cout << endl;

    }

}



int main()

{

    srand((unsigned)time(NULL));

    float **matrix = new float*[N];

    float **matrix2 = new float*[N];

    float **matrix3 = new float*[N];

    float **matrix4 = new float*[N];
    for (int i = 0; i<N; i++)

    {

        matrix[i] = new float[N];

        matrix2[i] = matrix[i];

         matrix3[i] = matrix[i];

         matrix4[i] = matrix[i];

    }



    //cout << "我们生成了初始随机矩阵" << endl;

    for (int i = 0; i < N; i++)

    {

        for (int j = 0; j < N; j++)

        {

            matrix[i][j] = rand() % 100;

        }

    }

    //print(matrix);



    cout <<endl<<endl<<endl<<"串行高斯消去法" << endl;

    clock_t  clockBegin,clockEnd;

    clockBegin = clock(); //开始计时

    float **B = normal_gauss(matrix);

    clockEnd = clock();

    //print(matrix);

    cout << "总共耗时： " << clockEnd - clockBegin << "ms" << endl;



    cout <<endl<<endl<<endl<< "使用SSE并行的高斯消去法" << endl;

    clockBegin = clock(); //开始计时

    SSE_gauss(matrix2);

    clockEnd = clock();

    //print(matrix2);

    cout << "总共耗时： " << clockEnd - clockBegin << "ms" << endl;



    cout <<endl<<endl<<endl<< "使用NEON并行的高斯消去法" << endl;

    clockBegin = clock(); //开始计时

    Neon_gauss(matrix3);

    clockEnd = clock();

    //print(matrix2);

    cout << "总共耗时： " << clockEnd - clockBegin << "ms" << endl;


    cout <<endl<<endl<<endl<< "使用AVX并行的高斯消去法" << endl;

    clockBegin = clock(); //开始计时

    AVX(matrix4);

    clockEnd = clock();

    //print(matrix2);

    cout << "总共耗时： " << clockEnd - clockBegin << "ms" << endl;


    system("pause");

    return 0;

}
